! ----------------------------
!
! Asks which window the user pretends to use
!
! ----------------------------
subroutine initPGPlot()
  implicit none
  
  integer*4 ier  
  integer*4 pgbeg

  ier = pgbeg(0, '/XSERVE', 1, 1)
  call pgask(.TRUE.)
  if (ier.NE.1) stop

  return
end subroutine initPGPlot


! ----------------------------
!
! Adjusts view to the given range and draws the axis and title labels.
! xMin   - minimum value for the axle Ox
! xMax   - maximum value for the axle Ox
! xMin   - minimum value for the axle Oy
! xMax   - maximum value for the axle Oy
! xLabel - label of the axle Ox
! yLabel - label of the axle Oy
! title  - label of the graph title
! axis   - used to specify the type of axis 
!          (correspond to the variable axis in the method PGENV)
!
! ----------------------------
subroutine setGraphView(xMin, xMax, yMin, yMax, xLabel, yLabel, title, axis)
  implicit none
  
  real*4 xmin
  real*4 xmax
  real*4 ymin
  real*4 ymax
  integer axis
  character*(*) xLabel
  character*(*) yLabel
  character*(*) title
    
  call PGPAGE
  call PGVSTAND
  call PGWINDOW (xmin, xmax, ymin, ymax)
  if (axis == 2) then
    call PGBOX ('BCGNTS', 0.0, 0, 'BCGNTV', 0.0, 0)
  else if (axis == 10) then
    call PGBOX ('BCNSTL', 0.0, 0, 'BCNTV', 0.0, 0)
  else if (axis == 20) then
    call PGBOX ('BCNST', 0.0, 0, 'BCNTLV', 0.0, 0)
  else if (axis == 30) then
    call PGBOX ('BCNSTL', 0.0, 0, 'BCNTLV', 0.0, 0)
  else
    call PGBOX ('BCNST', 0.0, 0, 'BCNTV', 0.0, 0)
  end if
  call PGBBUF
  call PGMTEXT('T', 2.0, 0.5, 0.5, title)
  call PGMTEXT('B', 3.2, 0.5, 0.5, xLabel)
  if (axis == 20 .or. axis == 30) then
    call PGMTEXT('R', 3.0, 0.5, 0.5, yLabel)
  else
    call PGMTEXT('L', 3.0, 0.5, 0.5, yLabel)
  end if
  call PGEBUF
  
  return
end subroutine setGraphView


! ----------------------------
!
! Computes the apropriate graph range for the given values
! and draws the axis and title labels
! X - vector with the x coordinates
! Y - vector with the y coordinates
! m - size of both vectors X and Y
! xLabel - label of the axle Ox
! yLabel - label of the axle Oy
! title  - label of the graph title
! axis   - used to specify the type of axis 
!          (correspond to the variable axis in the method PGENV)
!
! ----------------------------
subroutine setGraphViewAdjusted(X, Y, m, xLabel, yLabel, title, axis)
  implicit none
  
  real*4 X(m)
  real*4 Y(m)
  real*4 xmin
  real*4 xmax
  real*4 ymin
  real*4 ymax
  integer axis
  integer*4 m
  integer*4 n
  character*(*) xLabel
  character*(*) yLabel
  character*(*) title

  if (m > 0) then
    xmin = X(1)
    xmax = X(1)
    ymin = Y(1)
    ymax = Y(1)

    do n = 1, m
      if( X(n) < xmin ) then 
        xmin = X(n) 
      end if
      if( X(n) > xmax ) then 
        xmax = X(n) 
      end if
      if( Y(n) < ymin ) then 
        ymin = Y(n)
      end if
      if( Y(n) > ymax ) then 
        ymax = Y(n)
      end if
    end do
		
    if (xmin == xmax) then
      xmin = xmin - 0.5E0
      xmax = xmax + 0.5E0
    end if

    if (ymin == ymax) then
      ymin = ymin - 0.5E0
      ymax = ymax + 0.5E0
    end if

    !if (xmin > 0.E0) then
    !  xmin = 1.E0 * xmin
    !else if (xmin < 0.E0) then
    !  xmin = 1.E0 * xmin
    !else
    !  xmin = -xmax / 20.E0
    !end if

    !if (xmax > 0.E0) then
    !  xmax = 1.E0 * xmax
    !else if (xmax < 0.E0) then
    !  xmax = 1.E0 * xmax
    !else
    !  xmax = -xmin / 20.E0
    !end if

    !if (ymin > 0.E0) then
    !  ymin = 1.E0 * ymin
    !else if (ymin < 0.E0) then
    !  ymin = 1.E0 * ymin
    !else
    !  ymin = -xmax / 20.E0
    !end if

    !if (ymax > 0.E0) then
    !  ymax = 1.E0 * ymax
    !else if (ymax < 0.E0) then
    !  ymax = 1.E0 * ymax
    !else
    !  ymax = -ymin / 20.E0
    !end if
    
    call PGPAGE
    call PGVSTAND
    call PGWINDOW (xmin, xmax, ymin, ymax)
    if (axis == 2) then
      call PGBOX ('BCGNTS', 0.0, 0, 'BCGNTV', 0.0, 0)
    else if (axis == 10) then
      call PGBOX ('BCNSTL', 0.0, 0, 'BCNTV', 0.0, 0)
    else if (axis == 20) then
      call PGBOX ('BCNST', 0.0, 0, 'BCNTLV', 0.0, 0)
    else if (axis == 30) then
      call PGBOX ('BCNSTL', 0.0, 0, 'BCNTLV', 0.0, 0)
    else
      call PGBOX ('BCNST', 0.0, 0, 'BCNTV', 0.0, 0)
    end if
    call PGBBUF
    call PGMTEXT('T', 2.0, 0.5, 0.5, title)
    call PGMTEXT('B', 3.2, 0.5, 0.5, xLabel)
    if (axis == 20 .or. axis == 30) then
      call PGMTEXT('R', 3.0, 0.5, 0.5, yLabel)
    else
      call PGMTEXT('L', 3.0, 0.5, 0.5, yLabel)
    end if
    call PGEBUF
  end if
  
  return
end subroutine setGraphViewAdjusted


! ----------------------------
!
! Outputs the points in the current view
! X - array with the x coordinates
! Y - array with the y coordinates
! m - size of the arrays
!
! ----------------------------
subroutine outPoints(X, Y, m, pointStyle)
  implicit none

  real*4 X(m)
  real*4 Y(m)
  integer*4 m
  integer pointStyle

  if(m > 0) then
    call PGPOINT(INT(m), X, Y, pointStyle)
  end if

  return
end subroutine outPoints


! ----------------------------
!
! Outputs the points in the current view connected by lines
! X - array with the x coordinates
! Y - array with the y coordinates
! m - size of the arrays
! lineStyle - specifies the type of line to draw
!
! ----------------------------
subroutine outLines(X, Y, m, lineStyle)
  implicit none

  real*4 X(m)
  real*4 Y(m)
  integer*4 m
  integer lineStyle

  if(m > 0) then
    call PGSLS(lineStyle)
    call PGLINE(INT(m), X, Y)
  end if

  return
end subroutine outLines


! ----------------------------
!
! Outputs the matrix with gray scales. The lighter corresponds to the minimum 
! value and the darker corresponds to the maximum value.
! MATRIX - matrix with the values to display
! xSize  - the size of the matrix dimension corresponding to the x axle
! ySize  - the size of the matrix dimension corresponding to the y axle
!
! ----------------------------
subroutine outGrayMatrix(MATRIX, xSize, ySize)
  implicit none

  integer*8 i
  integer*8 k
  integer*8 xSize
  integer*8 ySize
  real*4 minValue
  real*4 maxValue
  real*4 TR(6)
  real*4 MATRIX(xSize, ySize)

  ! Initiates the min and max values
  minValue = MATRIX(1, 1)
  maxValue = MATRIX(1, 1)

  do i = 1, ySize
    do k = 1, xSize
      if (MATRIX(k, i) > maxValue) then
        maxValue = MATRIX(k, i)
      else if (MATRIX(k, i) < minValue) then
        minValue = MATRIX(k, i)
      end if
    end do
  end do

  TR    = 0.e0
  TR(2) = 1.e0
  TR(6) = 1.e0
  
  call PGGRAY(MATRIX, INT(xSize), INT(ySize), &
              1, INT(xSize), 1, INT(ySize), maxValue, minValue, TR)
  
  return
end subroutine outGrayMatrix


! ----------------------------
!
! Returns the array index in x axle and
! the corresponding real value in the y axle
! Y  - array with original values
! m  - size of the arrays
! XR - array to put the resulting x values
! YR - array to put the resulting y values
!
! ----------------------------
subroutine getValues(Y, m, XR, YR)
  implicit none
  
  real*8 Y(m)
  real*4 XR(m) 
  real*4 YR(m)
  integer*8 m
  integer*8 n

  do n = 1, m
    XR(n) = REAL(n)
    YR(n) = REAL(Y(n))
  end do

  return
end subroutine getValues


! ----------------------------
!
! Returns the array index in x axle and the Logarithm in 10th base in the y axle
! Y  - array with original values
! m  - size of the arrays
! XR - array to put the resulting x values
! YR - array to put the resulting y values
!
! ----------------------------
subroutine getLog10Values(Y, m, XR, YR)
  implicit none
  
  real*8 Y(m)
  real*4 XR(m) 
  real*4 YR(m)
  integer*8 m
  integer*8 n
  
  do n = 1, m
    XR(n) = REAL(n)
    YR(n) = REAL(DLOG10(DABS(Y(n))))
  end do
  
  return
end subroutine getLog10Values


! ----------------------------
!
! Returns the array index in x axle and the Natural Logarithm in the y axle
! Y  - array with original values
! m  - size of the arrays
! XR - array to put the resulting x values
! YR - array to put the resulting y values
!
! ----------------------------
subroutine getLogValues(Y, m, XR, YR)
  implicit none
  
  real*8 Y(m)
  real*4 XR(m) 
  real*4 YR(m)
  integer*8 m
  integer*8 n
  
  do n = 1, m
    XR(n) = REAL(n)
    YR(n) = REAL(DLOG(DABS(Y(n))))
  end do
  
  return
end subroutine getLogValues


! ----------------------------
!
! Returns the real part in x axle and the imaginary in y axle.
! U  - array of complex values
! m  - size of the arrays
! XR - array to put the resulting x values
! YR - array to put the resulting y values
!
! ----------------------------
subroutine getValuesFromComplex(U, m, XR, YR)
  implicit none
  
  real*4 XR(m)
  real*4 YR(m)
  integer*8 m
  complex*16 U(m)

  XR = REAL(DREAL(U))
  YR = REAL(DIMAG(U))
  
  return
end subroutine getValuesFromComplex


! ----------------------------
!
! Computes the slope and interception in y of the linear regression result.
! Only considers values between xMin and xMax, including.
! X - array with the x values
! Y - array with the y values
! xMin - minimum index of the points to consider
! xMax - maximum index of the points to consider
! m    - size of the arrays with the values
! slope      - returns with the computed slope
! yIntercept - returns with the interception of the rect with y 
!
! ----------------------------
subroutine getLinearRegressionData(X, Y, xMin, xMax, m, slope, yIntercept)
  implicit none
    
  real*4 X(m)
  real*4 Y(m)
  real*4 sxy
  real*4 sxx
  real*4 slope
  real*4 xAverage
  real*4 yAverage
  real*4 yIntercept
  integer*4 m
  integer*4 i
  integer*4 xMin
  integer*4 xMax
  
  sxy = 0.e0
  sxx = 0.e0
  xAverage = 0.e0
  yAverage = 0.e0
  
  do i = xMin, xMax
    xAverage = xAverage + X(i)
    yAverage = yAverage + Y(i)
  end do
  
  xAverage = xAverage / (xMax - xMin + 1)
  yAverage = yAverage / (xMax - xMin + 1)
  
  do i = xMin, xMax
    sxy = sxy + (X(i) - xAverage) * (Y(i) - yAverage)
    sxx = sxx + (X(i) - xAverage)**2
  end do
  
  slope = sxy / sxx
  yIntercept = yAverage - slope * xAverage
  
  return
end subroutine getLinearRegressionData


! ----------------------------
!
! Computes two points that define the rect resulted from a linear regression.
! Only considers values between xMin and xMax, including.
! X - array with the x values
! Y - array with the y values
! xMin - minimum index of the points to consider
! xMax - maximum index of the points to consider
! m    - size of the arrays with the values
! XLR  - array returns with the x values of the points computed
! YLR  - array returns with the y values of the points computed 
!
! ----------------------------
subroutine getLinearRegressionPoints(X, Y, xMin, xMax, m, XLR, YLR)
  implicit none
  
  real*4 X(m)
  real*4 Y(m)
  real*4 XLR(2)
  real*4 YLR(2)
  real*4 slope
  real*4 yIntercept
  integer*4 m
  integer*4 xMin
  integer*4 xMax
  
  call getLinearRegressionData(X, Y, xMin, xMax, m, slope, yIntercept)
  
  XLR(1) = X(xMin)
  XLR(2) = X(xMax)
  YLR(1) = slope * XLR(1) + yIntercept
  YLR(2) = slope * XLR(2) + yIntercept
  
  return
end subroutine getLinearRegressionPoints
